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Modeling near wall effects in second 
moment closures by elliptic relaxation 

By D. Laurence 1 and P. Durbin 2 


The elliptic relaxation model of Durbin (1993) for modeling near- wall turbulence 
using second moment closures (SMC) is compared to DNS data for a channel flow 
at Re t = 395. The agreement for second order statistics and even the terms in their 
balance equation is quite satisfactory, confirming that very little viscous effects (via 
Kolmogoroff scales) need to be added to the high Reynolds versions of SMC for 
near-wall-turbulence. The essential near- wall feature is thus the kinematic blocking 
effect that a solid wall exerts on the turbulence through the fluctuating pressure, 
which is best modeled by an elliptic operator. Above the transition layer, the effect 
of the original elliptic operator decays rapidly, and it is suggested that the log- 
layer is better reproduced by adding a non-homogeneous reduction of the return to 
isotropy, the gradient of the turbulent length scale being used as a measure of the 
inhomogeneity of the log-layer. The elliptic operator was quite easily applied to the 
non-linear Craft & Launder pressure-strain model yielding an improved distinction 
between the spanwise and wall normal stresses, although at higher Reynolds number 
(Re) and away from the wall, the streamwise component is severely underpredicted, 
as well as the transition in the mean velocity from the log to the wake profiles. In 
this area a significant change of behavior was observed in the DNS pressure-strain 
term, entirely ignored in the models. 


1. Introduction 

Second moment closures have the ability to account exactly for turbulence “pro- 
duction” terms due to shear, rotation and stratification, and to provide a better 
description than eddy viscosity models of the Reynolds stresses — a corner stone for 
complex flows involving heat transfer, two-phase flows, or combustion. However 
they are mainly used by industry in their high Re form, since near-wall models are 
both unsatisfactory and rather difficult to solve numerically. 

Since the publication of the budgets of the Reynolds stresses in a channel flow 
by Mansour, Kim & Moin (1988), a variety of near-wall second moment closures 
have been proposed (see review of 9 models by So et al., 1990). Some of them were 
more or less successful, but they are very seldom used outside low Reynolds number 
channel flows. Most of them use damping functions to force homogeneous models 
to comply with near wall turbulence features. A sound general principle is to avoid 
explicit use of the distance to the wall; however, this tends to render the models 
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rather difficult to converge numerically since the damping functions then depend 
on parameters of the unknown solution (such as the turbulent Reynolds number 
k 2 /ev). Launder &: Tselepidakis (1991), for instance, did a careful term by term fit 
to each component of the Reynolds stress balance obtained by DNS, but the overall 
model proved somewhat unstable especially when gravitational damping was added 
(Laurence, 1993). 

A common feature of these models is that all the terms (except diffusion) are 
related algebraically to the local values of the solution. Such local representation 
is in contradiction to the very large structures (hundreds of wall units) interacting 
with the wall and inhomogeneity of the velocity profile. 

In contrast to this previous approach, Durbin (1993) wets able to reproduce quite 
satisfactorily the features of near- wall flows by combining the very simple 4 IP’ ho- 
mogeneous second moment closure (Launder, Reece & Rodi, 1975) with a nonlocal 
(elliptic) approach representing the wall blocking effect on the large eddies. The 
present study was aimed at a closer comparison of the elliptic operator with DNS 
budgets and an analysis of what could be gained by combining it with a more 
sophisticated second moment closure (SMC). 

The Craft- Launder (1991) cubic SMC, either in free flows or in combination 
with wall functions or a low Re two-equation model, was shown to give better 
predictions than the IP model in a variety of flows (round and plane jets, impinging 
jets, tube bundles, swirling jets). A characteristic of near wall flows (also present in 
shear flows) is the strong reduction of the normal stress compared to the spanwise 
and longitudinal stresses. The Launder-Craft cubic model reproduces this effect 
in free flows and to some extent in the log-layer of a channel flow (Launder and 
Tselepidakis, 1991). Lee, Kim & Moin (1990) showed that many features of wall 
flows are also present in high shear homogeneous flows, though not accounting 
totally for the very high anisotropies in a near wall flow. Thus, our project was 
motivated by the idea that by combining the cubic model with the elliptic operator, 
to correctly acknowledge what is due to the high shear and what is due to the wall 
blocking effect, an improved model would result. 

2. Elliptic relaxation 

Following the procedure developed by Durbin (1993) the Reynolds-stress trans- 
port equation is written as: 

DtUiUj = Pij + Pij - + r «i + t / V 2 UjU J (1) 

Pij = -u t u k d k U] - u 3 u k d k Ui 

£ 

Pij = -uidjp - Ujdip - Eij + ttTuJ— (2) 

Tij = -d k u k UiUj 

The term pij differs from the usual pressure-strain (f> tJ since it includes the mis- 
alignment of the dissipation tensor and the Reynolds stress tensor: 
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Pij 



( 3 ) 


This unclosed term, called hereafter ‘relaxed pressure-strain’, is obtained by solv- 
ing an elliptic equation: 


£ 2 V2 £ll = 

k k 




( 4 ) 


For homogeneous turbulence p, j in Eq. 4 reduces to for which any standard 
redistribution model <p\j can be used: 


Pi j = 4>ij + dev(uiu j )^ 


( 5 ) 


where ‘dev’ is the deviatoric operator: 

dev(ujUj) = uiu~j — uWk £>ij/ 3 

The simple ‘IP’ model uses the Rotta return to isotropy and the ‘isotropization 
of production’; i.e., the slow and rapid parts are modeled as: 

= (p ij, slow 4" (pij, rapid 

F ( 6 ) 

= -C\ dev (uiUj)- — C% dev (Pjj) 

Durbin(1993) applied elliptic relaxation to the IP model with the following modifi- 
cations which define what is called hereafter the ‘R-linear model’: 

Pij = -(Ci - 1) deV ^ > -— - C 2 dev (Pij) (7) 


where the time scale is defined as: 


( 8 ) 


This time scale is also used in the dissipation equation in place of e/k, preventing 
a singularity at the wall. The length scale L appearing in (4) is also prevented from 
going to zero at the wall using again a Kolmogoroff scale as a lower bound: 


L = Cl max 



( 9 ) 


Furthermore Eq. 4 is solved numerically by introducing an intermediate variable 
fij = pij/k, and boundary conditions at the wall are imposed on the coupled 
u,Uj — fij equations. 

Last, the Daly-Harlow expression for the turbulent diffusion was used to model 
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y+ 

FIGURE 1, Reynolds Stresses, R-cubic model; □ , u\; o , a ? a*. Ujte2< 



FIGURE 2. Reynolds Stresses, R-linear model. Symbols as in fig. 1. 

3. Cubic order pressure-strain model 

Defining the anisotropy tensor and its invariants as: 

a \j — dev(ujUj)/k, A 2 — aijaij, A 3 = Qijajkdki, A = 1 — 9 / 8 (A 2 — A 3 ), (11) 
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FIGURE 3. Balance terms of lines, R-cubic model; symbols, DNS; , 

□ i p 22 • a , total diffusion; o , dissipation. 

the cubic pressure-strain model of Craft & Launder (1991) is written as: 

(praptd _ _ o.6d ev(Pij) + O.ZaijPkk 

- 0.2 +dkUi)~ -j^-( uJukdiUj + uJukdiUi)) (12) 

— r(Ai(Pij — Dij) + 3a m ,a„j(P mn — D m »*)) 

where: D tJ = uiutdjUk — ujuidiUk and 

4>ij° w — “ {Ci + l)a ( >£ (13) 

- C[Cidev(a ik akj)e 

Thus the R-cubic model is defined with the following expression for p on the RHS 
of (4): 

p% = -CM, + C;dev(a, *«*,))- + <t>T d ( 14 ) 

With the values: 

Ci = 3.1[Amin(A2,0.6)] 1 / 2 , C[ = 1.2, r = 0.6 

4. Low-Reynolds number channel flow 

Figs. 1 & 2 show the Reynolds stresses compared to the DNS data at Re=395 
(unpublished CTR simulation), as obtained by the relaxed cubic model and the 
relaxed IP model, respectively. The Uj component is slightly too close to the u\ 
component with the IP model. Both models predict a too steep decrease of u\ away 
from the wall. 
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Figure 4. Balance terms of u \ , R-linear model; captions as in Fig. 3. 


The budgets of the normal stress u\ are given in Figs. 3 and 4 for the R-cubic and 
R-linear models respectively. The dotted lines are the source term p f ■ in Eq. 4; the 
dashed line is the solution to Eq. 4 and provides the relaxed pressure-strain that 
enters the solutions shown in Figs. 2 and 3. Note that by y+ » 80, pij has relaxed 
to p^j. Figs. 3 and 4 are plotted on the same scale to show that the relaxed pressure- 
strain terms axe nearly identical even though the maximum of ^22 at y + = 20 in the 
homogeneous IP model is well out of range (0.11 in Fig. 4; i.e., 3 times that of the 
cubic model). This demonstrates that the solution to the relaxation equation has 
a large contribution coming through the boundary conditions. This blocking effect 
gives adequate near- wall behavior despite serious inaccuracies of the homogeneous 
model on the r.h.s of Eq. 4. The agreement with the balance terms obtained from 
DNS is satisfactory. For this comparison the DNS data were processed as in Eq. 1, 
in which Tij and uV 2 iTtU] are combined as total diffusion. Note that differences 
between the P 22 term and the corresponding DNS data are compensating for visible 
defects in the diffusion model, Eq, 10. Since the u\ profile was seen to be quite 
accurately predicted, the Daly Harlow model of turbulent diffusion needs to be 
revisited. 

Previous experience _with the standard Launder-Tselepidakis model exhibited dif- 
ficulties in solving the u\ balance because the normal stress component goes to zero 
as y 4 at the wall, while its balance terms remain large. In the form of Eq. 1, pressure 
strain now balances most of the diffusion while the remaining dissipation is a small, 
numerically stabilizing term. As — ► 0, the molecular diffusion balances the re- 

laxed pressure- strain P 22 = k f 22 , both of them going to zero, which is precisely the 
boundary condition imposed on the coupled system u\ — f 2 2 . 

A slightly less satisfactory agreement is obtained for the budget of uiu 2 shown 
in Fig. 5. The small overestimation of the normal stress seen in Fig. 1 results in 
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FIGURE 5. Balance terms of \ 1 \U 2 \ lines, R-cubic model; symbols, DNS.; o , 

production; , Pi 2 ‘, 0 » Pi 2 ; A > total diffusion; o , dissipation. 



FIGURE 6. Balance terms of u \ , R-cubic model; captions as Fig. 5. 

a visible overprediction of the production in the near wall layer, compensated by 
a similar overprediction of the pressure-strain. Since the extent of agreement for 
the latter is quite different depending on which component is examined, further 
improvement could only be obtained by a tensorial correction whereas only global 

coefficient tuning was undertaken here. 

Fig. 6 shows the balance of the streamwise stress u\. The production term be- 
comes dominant and the pressure-strain is now a small part of the budget. The 
dissipation is quite well predicted. The e model equation is only changed from its 
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standard high Re form by the use of T in place of k/e (Durbin, 1993). The pro- 
duction of dissipation needs to be slightly increased very near the wall. Hanjalic &; 
Launder (1976) introduced an extra production term, proportional to the second 
derivative of the mean velocity; but Rodi &: Mansour (1990) showed that such a 
term was too strong below y + = 10. In the R-linear model production is enhanced 
by modifying the production constant c ci as: 


z e x ~ C c\ + a i 


with a] 0.1. The dissipation equation is then: 


D t e = 


c' ei \Pkk ~C €2 £ 

T 


+ dk 





(16) 


A dependence on A 2 was added in the present R-cubic model to make sure the 
modification would have no effect on free shear flows, since A 2 only attains values 
near unity in the buffer layer where the turbulence becomes highly anisotropic: 


4, = c £l (l + aiming, l) 2 y) 


(17) 


The kink in £n near y+ = 10 is exaggerated but again compensates for a defect 

in the diffusion model. Note that in (16) no damping function is needed before c e2 
2 

since has been replaced by 4p which is finite. 

The distribution of the various structure parameters used throughout the paper 
are given in Fig. 7. The A 2 parameter becomes very large below y+ = 30, character- 
izing nearly 2-D turbulence and providing a means of isolating the transition layer, 
whereas the ratio of production over dissipation reaches values of about 1.5, which 
can also be found in free shear flows. The Rotta constant C\ from the Launder- 
Craft model goes to zero at the wall, indicating that one might need a smaller 
elliptic correction than the linear IP model. The difference c e2 — c ffl which can be 
related to the von Karman constant is seen to be fairly constant in the log-layer. 

We consider in Fig. 8 the split of <f> 22 into the slow part and the rapid part modeled 
by Eq. 12. The slow part (to which the quadratic component ($2) makes a small 
contribution) is dominating the rapid part and is acting ‘naturally’ to reduce the 
isotropy. The rapid part is determined largely by the linear term (the first term on 
the r.h.s. of Eq. 12, which is referred to as rj in figure 8) since the second and third 
terms (r2 and 7*3) seem to cancel each other. However, the model does not reduce 
to the linear contributions to the slow and rapid part ( s\ and ri) (which is just the 
IP model) because r3 takes a positive (‘natural’) sign in the spanwise stress budget. 
Also, in the shear stress budget r 2 is now acting ‘anti-naturally’ and 7*3 ‘naturally’, 
as can be seen in Fig. 9. 

Fig. 10 shows the velocity profile from the Re r = 395 DNS and the Comte-Bellot 
(1965) experiment at Re r = 2420 (Re = 57,000). The latter when matched with 
the standard log-law, {/+ = l/Klog(y+) + C, gives an additive constant of C=7, 
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y+ 

FIGURE 7. Structure parameters; see Eqs. 11 and 14. a , ce 2 — ct\; o, A; a , A 2 ; 
«, Cli V, P/e; *, dl/dy . 



y+ 

Figure 8. Split of (f> 22 in the cubic model, o , ^r, □ , s 2 ; o, rj; A, r 2 ; <3, r 3 ; 
, total. 

which is unusually high but will serve here to emphasize the effect of Ct in Eq. 9. 
Models are usually calibrated to yield C=5.5 and /c=.41 (dot-dashed line). The 
constant C is known to have a Reynolds dependence and is related to the Van 
Driest damping factor A+\ i.e., it depends on the rate at which the shear stress 
increases in the transition in the region 10 < y+ < 30, effectively accelerating the 
mean flow. This also introduces a pressure-gradient dependence. 

With the R-cubic model described up to now, a best fit with the log-profile was 
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o.o 


50.0 100.0 150.0 200.0 


FIGURE 9. Split of 4> 12 in the cubic model. Symbols as in fig. 8. 


obtained with Cl — 0.25. The homogeneous version of the cubic pressure-strain 
model yields an overestimation of the slope of the velocity profile. Thus Launder- 
Tselepidakis included a Gibson- Launder type of wall reflection (Eq. 24) to reduce 
the pressure-strain return to isotropy. Alternatively, they introduced an effective 
velocity gradient in the pressure strain model defined as: 


l 1 / 2 tj — = 7 — 

VUi eff = VUi + c e// /(V/ • V)V(t/,) , l = - — 


In the present 1-D problem, this is equivalent to: 

dU eff =dU eff'dj&U 

dy dy + ° dy dy 2 

Or, since <^ J>r is linear in the mean shear: 


4> e !f = fa T + l — ^ 
9%i ' r 9,1 ’ r ^ c dy dy 


(18) 


(19) 


( 20 ) 


In the log-layer <j> behaves as y~ l so the effect is to reduce <f>ij by a factor which can 
be estimated using standard log-law assumptions as about 40%. 

Note that this non-local effect could be incorporated in the relaxation operator 
by replacing kL 2 V 2 ^- in (4) by: 


V-(I 2 V^) = I 2 V 2 ^ + 2IV(L)v£2 (21) 

In the log region, if we assume py = <f>i j, and neglect the effect of variations of A, 
the first term on the RHS has the same sign as pij thus increasing the return to 
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y+ 


FIGURE 10. Velocity profiles, R-cubic model. , c{ ff = 0.12, — 0.1, cl 

= 0.15; , clefF = 0.12, c2eff =0.1, cl = 0.22; , U+ = 1/k log(y+) + C; 

o , DNS; □ , Comte- Bello t 

isotropy (observed if Fig. 4 is plotted for y+ > 100) while the second term actually 
reduces the return to isotropy, which is the effect sought by using wall echo terms. 

This interesting idea, which avoids any explicit reference to the distance to the 
wall, was for the time being retained in a simpler version as follows: 



The following combination: 

c\ fS = .12, c e 2 ff = .1, C L = .15 

yielded the more satisfactory agreement, increasing U in the transition layer and 
slightly decreasing the slope of U further away from the wall in better agreement 
with the log-law as shown by the dashed curve in Fig. 10. The solid line obtained 
with c\ ff = .12, c\^ = .1, Cl = .22 shows that in the present form, the relaxation 
effect is limited to the transition layer, whereas in the R-linear model used without 
Eqs. 22, 23, it also strongly affects the log-layer (Fig. 11). 

The final form of the model at this stage and for which the results were shown 
in Figs. 1-10 is thus the standard Craft-Launder model, with the elliptic relaxation 
and the modifications in Eq. 22 and the values: 

c £l c £j c M a e Ok Cl C n ai c\ ff c\ ss 

1. f(A,A 2 ) .23 1.3 1. .15 80. .2 .12 .1 
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V+ 

FlGURE 11. Mean velocity profiles, R-linear model. For captions see Fig. 10 

while for the R-linear model used here the constants of Durbin (1993) were used: 

Cp &k Cl Cjj CLi C\ C 2 

1.44 1.9 .23 1.65 1.2 .2 80. .1 1.22 0.6 

Note that the pressure strain constant C\ is lower than the standard value; this 
reduction could be avoided by using the gradient of the length scale as an inhomo- 
geneity indicator. Also the following Gibson- Launder formulation, 

<f>Tj = <t>kinkrn6ij - 3/2 4>ikn k rij - 3/2<£jjtn*n;, (24) 

was not used with the Craft-Launder model since u\ is sufficiently suppressed by 
the rapid term of the cubic model. 

At Re=395, Fig. 12 shows that the R-linear and R-cubic predictions are almost 
undistinguishable. The fact that all models known to the authors seem to under- 
predict the increase in velocity in the central part of the channel (the wake region) 
where they recover their homogeneous form is somewhat puzzling. The centerline 
velocity (and more importantly, the skin friction in boundary layers) seems to be 
recovered only at the expense of predicting a somewhat lower von Karman constant. 
Considering the UiU 2 balance equation, neglecting diffusion and dissipation, 


0 = -n 2 2 -^+<f>i2, (25) 

one sees that with u\ constant, the magnitude of the velocity gradient is allowed 
to increase (relative to y -1 in the log layer) only if the pressure strain decreases 
less than y" 1 . This is indeed the behavior of exhibited by the DNS data (Fig. 13) 
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Figure 12. Velocity at Re=395: , R-cubic; , R-linear. For captions 

see Fig. 10. 



y+ 

FIGURE 13. Budget of terms in the u 2 equation near channel centerline, captions 
as Fig. 4. 

and seems to be ignored by the models. In fact the relative increase of the pressure 
strain exceeds that of the production and is balanced by pressure diffusion. A 
similar behavior was found for all Reynolds-stress budgets. 

With the present gradient transport assumption, a zero value is predicted for 
the total diffusion (pressure + turbulent) since iTjUJ is linear, but non-zero pres- 
sure diffusion might be accounted for by the non-local formulation investigated by 
Demuren et al in the present volume. 
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FIGURE 14. Reynolds stresses for the channel flow at Re=57,000. o , ti|; o , 

□ , 1 * 3 ;^, u\u 2 . , R-cubic; , R-linear. 

Fig. 14 compares the Reynolds-stress distributions throughout the channel at high 
Reynolds number (Comte-Bellot experiment Re = 57,000). One major drawback 
is that u\ is severely underestimated. Experiments show an extension with Re of 
the plateau region of the stress maxima (Antonia et a/., 1992), whereas the models 
follow this trend more moderately. This is a generic problem since the standard 
high Re Gibson- Launder model with wall functions yields similar underestimations 
away from the wall. One sees, however, that the R-cubic model reproduces a better 
separation between normal and spanwise components than the R-linear model for 
which this separation is limited to the log layer. 

5* Other calculations 

The skin friction was computed for a zero pressure gradient, boundary layer. The 
R-cubic model overpredicts C/, which is related to the difficulty in predicting the 
wake region in the centerline of channel flow that was observed previously (Fig. 10). 

The R-cubic model was also tested on the flow over a backward facing step but 
without Eqs. 22-23 and yielded similar, if not less satisfactory, predictions in com- 
parison to the R-linear model used by Ko &; Durbin (1993). 

Conclusion 

The present study used the DNS results of a channel flow at Re r = 395 to confirm 
that homogeneous, second moment closures can be quite easily made to comply with 
near- wall turbulence characteristics by applying the elliptic relaxation procedure of 
Durbin (1993). Physically, it models the blocking effect that the wall imposes on 
the the fluctuating pressure, thus alleviating the need for Re dependent ‘damping 
functions’. 
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It was found that the standard relaxation model produces a reduction of ‘return 
to isotropy’ in the near wall layer (y+ < 80). This effect is essentially due to the 
boundary conditions. In this region, the elliptic relaxation is so strong that the 
switch from a linear to a cubic pressure strain model had a nearly unnoticeable 
effect on the budgets of the stresses. After imposing the elliptic relaxation these 
budgets compare very well with the DNS data. An urgently needed improvement 
concerns the Daly-Harlow turbulent diffusion term which was not studied here. 

Further away from the wall, it seems that the strong inhomogeneity of the log- 
layer has also a significant blocking effect, underestimated by the original elliptic 
relaxation combined with the simple IP second moment closure. Using the Craft & 
Launder model, the wall normal stress was better reproduced, but still insufficiently 
to avoid further inhomogeneity corrections. This is the reason that the Gibson 
Launder ‘wall echo’ model was still required at significant distances from the wall. 
In the latter, reference to the distance to the wall can be avoided by using the 
gradient of the turbulence length scale as suggested by Launder & Tselepidakis, 
and which could be included along with elliptic relaxation. 

Outside the log-layer, the DNS data show a significant change in the behavior of 
the pressure strain terms, which explains the increase of the velocity gradient, but is 
not reproduced by the models. The latter seem to compensate for this omission by 
overpredicting the slope in the log-layer. Near wall models are usually compared to 
DNS data at low Re, but for practical applications more attention should be given 
to higher Re flows, with the challenging feature that the profiles of the stresses show 
a strong Re dependence (i.e., they do not collapse on plots seeded in wall units). In 
particular the streamwise stress is severely underpredicted at high Re. 
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